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Robustness to mutations and noise has been shown to evolve through stabilizing selection for 
optimal phenotypes in model gene regulatory networks. The ability to evolve robust mutants is 
known to depend on the network architecture. How do the dynamical properties and state-space 
structures of networks with high and low robustness differ? Does selection operate on the global 

OO ' dynamical behavior of the networks? What kind of state-space structures are favored by selection? 

^P I We provide damage propagation analysis and an extensive statistical analysis of state spaces of these 

model networks to show that the change in their dynamical properties due to stabilizing selection for 
optimal phenotypes is minor. Most notably, the networks that are most robust to both mutations 
?— j I and noise are highly chaotic. Certain properties of chaotic networks, such as being able to produce 

large attractor basins, can be useful for maintaining a stable gene-expression pattern. Our findings 
indicate that conventional measures of stability, such as the damage-propagation rate, do not provide 
much information about robustness to mutations or noise in model gene regulatory networks. 
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I. INTRODUCTION 

_^ ■ The genetic architecture of biological organisms shows remarkable robustness against both structural and envi- 

^ ' ronmental perturbations [71, [ill, [iTI, |45|. For example, quantitative models indicate that the functionality of the 

T^ • Drosophila segment polarity gene network is extremely insensitive to variations in initial conditions [4]| and robust 

against architectural modifications [42|. (Also see Albert and Othmer [l|], Chaves et al. [3|-) Gene knock-out studies 

on yeast have shown that almost 40% of the genes on chromosome V have either negligible or no effects on the growth 

rate [33|. Certain cellular networks, such as the E. Coli chemotaxis network, are also known to be very robust to 

\jj \ variations in biochemical parameters [J, l7||. 

J^ • The genetic regulatory networks that control the developmental dynamics buffer perturbations and maintain a stable 

pg \ phenotype. That is why phenotypic variation within most species is quite small, despite the organisms being exposed 
l/~j , to a wide range of environmental and genetic perturbation s Usll. I t has been proposed that genetic robustness evolved 
through stabilizing selection for a phenotypic optimum [111 . l43l . |44| . Wagner [43| showed that this in fact can be true by 
modeling a developmental process within an evolutionary scenario, in which the genetic interaction sequence represents 
the development, and the stationary configuration of the gene network represents the phenotype. His results indicate 
t — \ that the genetic robustness of a population of model genetic regulatory networks can gradually increase through 
^D ■ stabilizing selection, in which deviations from the "optimal" stationary state (phenotype) are considered deleterious 
-- 0. 

In this paper, we focus on the effects of evolution of genetic robustness on the dynamics of gene regulatory networks 
in general. First, we examine the relationship between genetic robustness and the dynamical character of the networks. 
By dynamical character, we mean the stability of the expression states of the network against small perturbations, or 
noise. 

Models that are employed to study dynamics of gene networks, such as random Boolean networks (RBN) [2, l3, llSl, 
122 . |23| or some variants of random threshold networks (RTN) [26| have been known to undergo a phase transition at 
low connectivities [30, [M|' giving rise to change in dynamical behavior: on average, small perturbations will percolate 
through the network above the threshold connectivity (chaotic phase), whereas they stay confined to a part of the 
network below the threshold (ordered phase) . Intuitively, one might expect to find that robustness to mutations (which 
are permanent structural changes, not dynamic perturbations) is related to the dynamical behavior of the system, 
therefore, ordered gene regulatory networks should be genetically more robust than the chaotic ones. However, this 
is not necessarily true. Here, we show that the relation between the dynamical character of a genetic regulatory 
network and its mutational robustness can be quite the opposite. In fact, even earlier studies provide a clue on this 
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issue: for gene networks that have undergone selection (i.e., evolved), mutational robustness is known to increase with 
increasing connectivity [43|- On the other hand, chaoticity has been shown to increase with increasing connectivity 
in random gene regulatory networks 0, [3J|. These facts seem to contradict the intuitive interpretation of robustness 
since they suggest that chaotic networks can be mutationally robust. 

Although this inference is proven true in this paper, such an assessment cannot be based on the previous studies 
of the dynamics of random networks (3| . pa . [25l. [2a. [sol. [34|. as the evolved networks (with high mutational robustness 
and connectivity) mentioned above [43| have undergone selection. The selection process could potentially tune a 
network to exhibit a different dynamical character than its random ancestors. Therefore, the dynamical character of 
the evolved networks should be studied independently and then compared with their random counterparts. Here, we 
study the dynamics of the evolved networks numerically and show that selection for an optimal phenotype indeed 
has only a minor effect on their global dynamical behavior. This result indicates that the evolution of mutational 
robustness cannot be understood in terms of simple dynamical measures. 

We also provide statistics on robustness to noise, which is the ability of a network to reach its "optimal" steady state 
after a perturbation to the gene-expression trajectory. Computer simulations indicate that mutational robustness is 
correlated to robustness of the gene-expression trajectory to small perturbations (noise), at least for short trajectories. 
This result is supported by recent studies [id. l2l|. For perturbations of arbitrary magnitude, the basin size of the 
steady-state attractor provides a better measure. Our analysis shows that basin sizes of densely connected networks 
(which are highly chaotic) have a broad distribution, and therefore such networks can have very large attractor basins 
[3, 24, 28]. The yeast cell-cycle network has been reported to have similar properties [2y|. Although chaoticity is just 
a side-effect of high connectivity, and not directly selected for during evolution, it appears that this intrinsic property 
of chaotic networks can be useful in terms of robustness to noise. When all of these measures are taken into account, 
chaotic dynamics does not seem be an obstacle for the networks with high connectivity since they are more robust to 
both mutations and noise than the ones that are sparsely connected. 

The organization of the rest of this paper is as follows. We describe the model in SecHIland explain its implemen- 
tation in simulations in Sec. lIIII We give the results in Sec. HV] and discuss their implications in SeclVl 

II. MODEL 

We use the model introduced by Wagner [43|, which has also been used with some modifications by other researchers 
0, S [3fl|. Each individual is represented by a regulatory gene network consisting of N genes. The expression level 
of each gene, s^, can be either -1-1 or —1, meaning that the gene is expressed or not, respectively. The expression 
states change in time according to regulatory interactions between the genes. The time development of the expression 
states (i.e., the dynamical trajectory taken by the network) represents a developmental pathway. This deterministic, 
discrete-time dynamics of the development is given by a set of nonlinear difference equations, 

s,(t + 1) = I ^^"^ {^f=i w^jSj{t)j , J2f^^ w,jsj{t) ^ 

where sgn is the sign function and w^ is the strength of the influence of gene j on gene i. Nonzero elements Wij of 
the N X N matrix W are independent random numbers drawn from a gaussian distribution with zero mean and unit 
variance. The diagonal elements of W are allowed to be nonzero, corresponding to self-regulation. The mean number 
of nonzero elements in W is controlled by the connectivity density, c, which is the probability that any given Wij is 
nonzero. Thus, the mean degree of the network is {k) = cN. To clarify, t represents the developmental time, through 
which genetic interactions occur. It is different from the evolutionary time, T, which will be explained in the next 
section. 

The dynamics given by Eq. ^ can display a wide variety of features. For a specifled initial state s(0), the network 
reaches an attractor (either a fixed point or a Hmit cycle) after a transient period. Transient time, number of attractors, 
attractor periods, etc. can differ depending on the connectivity of the network, and from one realization of W to 
another. The fitness of an individual is defined by whether it can reach a developmental equilibrium, i.e., a fixed 
point, which is a fixed gene-expression pattern, s*, in a "reasonable" transient time. (It has been shown that selection 
for developmental stability is sufficient for evolution of mutational robustness [33|, i.e., deviations from s* do not have 
to be deleterious, as long as the gene-expression configuration reaches a fixed point. However, we shall adopt the 
criterion used by Wagner [43| for compatibility.) Further details of the model are explained in the next section. 



III. METHODS 
A. Before evolution: Generation of networks 

We studied populations of random networks with N = 10. We use these relatively small networks to be able to 
enumarate all network states exhaustively for a large set of realizations. In the simulations, each network was first 
assigned a random interaction matrix W and an initial state s(0). W was generated as follows. For each Wij , a random 
number uniformly distributed on [0, 1) was generated, and Wij was set to zero if the random number was greater than 
the connectivity density, c. Otherwise, Wij was assigned a random number drawn from a gaussian distribution with 
zero mean and unit variance. Then, each "gene" of the initial configuration, Si{0), was assigned either —1 or +1 at 
random, each with probability 1/2. 

After W and s(0) were created, the developmental dynamics were started, and the network's stability was evaluated. 
If the system reached a fixed point, s*, in 3N time steps, then it was considered viable and kept. Otherwise it was 
considered unstable, both W and s(0) were discarded, and the process was repeated until a viable network was 
generated. (A viable network can have several fixed points and/or limit cycles in addition to s*.) For each viable 
network, its fixed point, s*, was regarded as the "optimal" gene-expression state of the system. This is the only 
modification we made to the model used by Wagner ^3|: we accept any s* as long as it can be reached within 3A^ 
time steps from s(0), whereas Wagner [43| generated networks with preassigned random s(0) and s* [46|. 

B. Evolution 

In order to generate a collection of more robust networks, a mutation-selection process was simulated for each viable 
network as follows. First, a clan of A/" = 500 identical copies of each network was generated. For each member of the 
clan, a four-step process was performed for T — 400 generations: 

1. Recombination: Each pair of the N rows of consecutive matrices in the clan were swapped with probability 1/2. 
Since the networks were already shuffled in step 4 (see below), there was no need to pick random pairs. 

2. Mutation: Each nonzero Wij was replaced with probability l/{cN'^) by a new random number drawn from 
the same standard gaussian distribution. Thus, on average, one matrix element was changed per matrix per 
generation. 

3. Fitness evaluation: Each network was run starting from the original initial condition s(0). If the network reached 
a fixed point, s''', within t = 3N developmental time steps, then its fitness was calculated using 

f(st,s*)=exp(-H2(st,s*)/a,)), (2) 

where H(s^,s*) denotes the normalized Hamming distance between s''^ and s*, CTs denotes the inverse of the 
strength of selection, s* is the optimal gene-expression state, which is the final gene-expression state of the 
original network that "founded" the clan. We used CTs = 0.1. If the network could not reach a fixed point, it was 
assigned the minimum nonzero fitness value, exp(— l/o-g). 

4. Selection/Asexual Reproduction: The fitness of each network was normalized to the fitness value of the most 
fit network in the clan. Then a network was chosen at random and duplicated into the descendant clan with 
probability equal to its normalized fitness. This process was repeated until the size of the descendant clan 
reached J\f. Then the old clan was discarded, and the descendant clan was kept as the next generation. This 
process allows multiple copies (offspring) of the same network to appear in the descendant clan, while some 
networks may not be propagated to the next generation due to genetic drift. 

At the end of the T ~ 400 generation (evolutionary time) selection, any unstable networks were removed from the 
evolved clan. 

Some steps of the process above may be unnecessary, and in fact, the results do not depend strongly on model 
details [la, |43|- Nevertheless, the entire procedure of Wagner [43| was retained for compatibility. 

C. Assessment of Robustness 

The mutational robustness, i?^, of a network was assessed as follows. First, a nonzero Wij was picked at random and 
replaced by a new random number with the same standard gaussian distribution. Then, the developmental dynamics 
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FIG. 1: Probability density of mutational robustness (-R^) of 10000 sample networks before (filled bars) and after (empty bars) 
evolution with iV = 10 and (fc) = iVc = 10. Before: (7?^) = 0.63, a = 0.16. After: (i?^) = 0.92, a = 0.06. The evolved 
distribution was calculated by sampling one network from each of 10000 evolved clans. The mean indegree, (k), rather than 
the connectivity density, c, is the parameter that controls the behavior of the system [3, ll^l- See Sec. IIIII for other model 
parameters. 

were started, and it was checked whether the system reached the same stationary state, s*, within t = 3N time steps. 
This process was repeated 5000c times, starting from the original matrix. The robustness of the original network 
before evolution was defined as the fraction of singly- mutated networks that reached s*. For the evolved networks, we 
picked one sample network at random from the clan and used the same procedure to assess its mutational robustness. 

IV. RESULTS 

The stabilizing selection described above increases the robustness of a model population of gene networks against 
mutations [43|. Figure [T] shows that a population with a very large initial variation in robustness evolved an increased 
ability to absorb mutations after T = 400 generations of stabilizing selection. However, it is not very clear what kind 
of a reorganization in the state spaces of these networks occur during the evolution. We measured the changes in 
several system parameters to answer this question. 



A. Change in Dynamical Character 



It has been analytically shown that the RTNs undergo a phase transition from order to chaos with increasing (k) 
[34 |. Obviously, the state space of an RTN is finite. Therefore, the expression states have to display periodicity in 
development after at most 2^ steps. Thus, chaos here is not a long-term aperiodic behavior; rather it corresponds 
to a dynamical regime in which small perturbations percolate through the gene network [3|, l30|, l3]| . We quantify the 
dynamical character of a network by comparing the time development of two configurations, s(t) and s'(i), that differ 
by one gene: we measure the mean number of different genes in time step t+1, (dt+i), averaged over all possible _s(i) 
and s'(i) pairs with H(s(t), s'(i)) = 1. This is known as damage spreading or damage propagation|47| |l2l . [30 
The dynamical behavior of random networks can be quantified analytically usingdamage-spreading analysis [3C 
However, this analytical approach may not be applied to the evolved networks [43| as the selection process can tune 
a network to behave dynamically very differently than its ancestor. Therefore, we calculated (dt+i) numerically for 
both viable and evolved networks, averaging ensembles of 10000 networks exhaustively over all possible configuration 
pairs, s(t) and s'(t), for each (k). As seen in Fig. [21 the (dt+i) increases monotonically with increasing connectivity, 
indicating that highly connected networks are more chaotic on average. The evolved networks are slightly more 
ordered (on average) than their viable ancestors. These results indicate that the dynamical behavior of the evolved 
networks is not much different than that of the viable networks from which they are descended. 
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FIG. 2: Mean damage spreading, (rft+i), after a one-bit perturbation before (viable) and after evolution (evolved) for networks 
with A^ = 10, measured exhaustively for all possible state pairs, averaged over 10000 samples each. The error bars represent 
one standard deviation (not standard error). Evolved networks (triangles) are slightly more ordered compared to their viable 
ancestors. (The differences are statistically significant.) (dt+i) is always larger than unity due to our specific update rule. See 
the discussion in the text for details. The lines connecting the symbols are guides to the eye. (Probability distributions for 
(dt+i) for certain values of {k) are given in supplementary Figs. S3 and S4.) 

B. Change in the State-Space Structure 

We also analyzed the state spaces of random, viable, and evolved networks [33|- Figures [4] and [5] show these statistics 
for networks with A^ = 10 and connectivities (k) =1, 5, and 10. Due to the up-down symmetry of the system, the 
state space is divided into two parts, where the dynamics on one side is the mirror image of the other. Therefore, 
the basin size of a fixed-point attractor cannot exceed half the size of the state space, Q/2 = 2^~^. Limit cycles, 
however, can cross the symmetry plane that divides the state space (Fig. [Sljb)). Therefore, their basins can contain 
up to J7 = 2^ states. 

Typically, sparsely connected gene networks have many attractors (fixed points and limit cycles) and, consequently, 
smaller basin sizes on average, as depicted in Fig.jSja). Increasing connectivity brings state spaces with fewer attractors 
(Figs. UK a), (b), and (c)), longer attractor (limit cycle) periods (Figs.|Hd), (e), and (f)), and broadly distributed basin 
sizes (Figs,[4ljg), (h), and (i)). Densely connected networks can contain basins that occupy a large portion, occasionally 
even all, of the state space. The basins of such networks tend to have broadly distributed branch lengths as seen in 
Fig. [Sljb), and their states also tend to have fewer precursors (Figures [UJa) , (b), and (c)). (All states, {s(i)}, that go 
to state s''' at time step t+1 are precursors of s'''.) Therefore, mean transient time (total number of steps from a state 
to the attractor) on such networks are typically larger than on networks with lower connectivity (Figs.[5l[d), (e), and 

(f))- 

These changes in the state-space characteristics are consistent with the damage-spreading measurements. The dis- 
tributions for the evolved networks are shifted toward the distributions of more ordered networks of lower connectivity. 
This means evolved networks are slightly more ordered, as the damage-spreading measurements shown in Fig. ^h) 
suggest. For example, the distribution of the number of attractors (Figs.IDJa), (b) and (c)) for the networks with low 
connectivity ((fc) = 1) has a longer tail compared to the networks with (k) = 5 and 10. Similarly, the distributions for 
the evolved networks have slightly longer tails compared to those of viables of the same connectivity. The same effect 
can be seen in the attractor-period (Figs.SJd), (e), and (f)), precursor (Figs. [DJa), (b), and (c)), and transient-time 
distributions (Figs.[5jd), (e), and (f)), as well. The most significant changes are seen in the transient-time distribu- 
tions, indicating that basins with relatively shorter branches are preferred by selection. The basin-size distributions 
of viable and evolved networks, however, do not display much difference as they virtually overlap. 

C. Change in Gene Expression Trajectory and the State-Space Structure Nearby 



As seen in Fig.[6lja), networks of all connectivities have similar values of mutational robustness, i?^, before selection. 
After selection, however, the networks with higher connectivity reach a much greater mean mutational robustness 




FIG. 3: Two typical state-space structures shown as directed graphs. Each node represents a gene expression state, and the 
links represent developmental transitions from f to t + 1. Arrows indicate the direction of flow. Each state drains into an 
attractor (pentagons) through transient states, (a) State space of a sample network with N = 8 and (fc) = 1. There are eight 
basins, each having a flxed-point attractor. The principal basin (at the top) contains the initial (diamond shaped node) and 
flnal states. Note that the basins are quite symmetric, and transients are very short as the system behaves less chaotically 
when the connectivity of the network is low. (b) Basin of a sample network with A'' = 8 and (k) = 8. The principal basin 
occupies a large portion of the state space. The basin on the left crosses the symmetry plane with a 2-cycle (note the symmetry 
of the branches). The other basins have mirror images on the other side of the state space (not shown). The principal basin 
is signiflcantly larger than the others. The high connectivity makes the network behave more chaotic, creating basins with 
broadly distributed sizes and branch lengths, and significantly longer transients. See text for details. Generated using Graphviz 



[43|- Sparsely connected networks do not show much improvement because the low connectivity leaves little room 
for optimization. Another significant change is seen in the transient time, t (Figs. [DJg), (h), and (i) and[6l[b)) |43l |. 
For viable networks, the transient time increases with the connectivity. This is not surprising since state spaces of 
densely connected networks typically have basins with longer branches (Fig. [3ljb)). However, selection brings the 
mean transient time down to around 2, independent of the connectivity of the network. Low transient time is one of 
the properties of highly robust networks [10| . 

Mutations are not the only kind of perturbations that genetic regulatory networks experience. All genetic systems 
are also exposed to noise created by both internal an external sources. We measure robustness to noise - dynamical 
robustness - using two parameters. The first one is the dynamical robustness, i^, of a gene expression state, s, to 
small perturbations (random one-bit fiips) , which is simply the fraction of nearest neighbors in Hamming distance of 
s that lie in the same basin as s [10||. In other words, v is the probability that a random fiip of a gene on s yields a 
state that drains into the same attractor as s. We also employ the mean robustness of a set of states, for instance, 
robustness of the gene expression trajectory, vt, which is the average of ly over all states in the trajectory (including 
s(0) and s*). Similarly, the robustness of the principal basin, t'B, is the average of v over all states in the basin, and 
the robustness of the entire state space, vg, is the average of v over all possible states of the network. (We usually 
omit the term "dynamical" when we talk about robustness of a state or a set of states since the mutational robustness 
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FIG. 4: Probability distributions for the attractor-count (first ro-w), attractor-period (second ro-w), and basin-size (all basins, 
third ro-w) distributions for A'' = 10 and (k) = 1, 5, and 10 (columns 1, 2, and 3, respectively). Each curve represents an average 
over 20,000 realizations for the "random" net-works, and 10,000 realizations for the "viable" and "evolved" ones. For the evolved 
net-works, we did not use clan averages to avoid a bias. Instead, we picked one sample net-work from each evolved clan. Error 
bars -were calculated by grouping the data: Each data set -was divided into groups of 1000 samples and the average distribution 
for each group -was calculated. Then, the error calculations -were performed on the ne-w set of averaged distributions. The data 
plotted on log-log scale -were histogrammed using exponential bins (0, 1, 2-3, 4-7,...) to reduce the noise. (See supplementary 
Fig. SI for histograms -with linear bins.) The basin size does not include the size (period) of the attractor. The evolved 
attractor-count and attractor-period distributions are shifted to-ward those of viable net-works -with lo-wer (k) , indicating that 
evolved net-works display a slightly more ordered character. The basin-size distributions for the evolved net-works seem to 
overlap -with the ones for the viable net-works except for the largest basins, -which become less probable after selection. The 
lines connecting the symbols are guides to the eye. 



of a state is not defined.) 

The second parameter -we use to estimate dynamical robustness is the normaHzed size of the principal basin, B, 
■which is the basin containing both s(0) and s*. The size of a basin is the total number of states it contains. Therefore, 
B is the fraction of all possible perturbations to the trajectory that leave s* unchanged, i.e., the probability that 
flipping an arbitrary number of genes at random in a state on the trajectory does not change the attractor that the 
network settles into. 

Clearly, the parameters -we employ to measure robustness to noise are quite different than conventional measures of 
stability, such as the damage-spreading rate. This is because -we are interested in the endpoints of the gene expression 
trajectory, but not in the exact path taken. Therefore, the relationship bet-ween the damage spreading and robustness 
to noise is not trivial. In fact, it is quite counterintuitive as explained belo-w. 

As sho-wn in Figs. IHg), (h) and (i), densely connected net-works have a greater number of larger basins as their 
basin-size distributions essentially follo-w a power law for about two decades. The effect of the connectivity on B 
is similar, as shown in Fig. [6l[c). Networks with low connectivity have small principal basins both before and after 
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FIG. 5: Probability distributions for precursor (indegree of the state-space graph) distribution (first row), transient time for 
all states (second row), and transient time for the initial state (third row) for A'' = 10, (fc) = 1, 5, and 10 (columns 1, 2, and 
3, respectively). The second row shows the distributions for the lengths of the transients starting from all possible states. The 
third row, however, shows the distribution for r, the length of the transient starting from s(0), the original initial state. Again, 
changes in the distributions indicate that evolved networks display a slightly more ordered character. The data plotted on 
log-log scale were histogrammed using exponential bins (0, 1, 2-3, 4-7,...) and the other plots were histogrammed using linear 
(0,1,2...) bins. See the caption of Fig. |4]for the details of the data analysis. The lines connecting the symbols are guides to the 
eye. 



selection. More densely connected networks have larger B before selection. We see a small increase in networks with 
(fc) = 5, while fully connected networks ((fc) = 10) increase their (B) significantly after evolution. However, compared 
to the change in mutational robustness, the effect of selection on the size of the principal basin is modest for highly 
connected networks. 

Although greater connectivity increases the damage spreading (sensitivity to small perturbations), the dynamical 
robustness of the trajectory, i^t, as well, increases with the connectivity (Fig. [6ljd)). Networks with (fc) = 1 have a 
small (i^t) both before and after selection. The V'^ curve for the viables has a maximum around (fc) = 2.5, which 
also shows no improvement after selection. This indicates that the networks with (/c) « 2.5 are intrinsically robust 
to noise, but not evolvable. Viable networks with higher connectivity have monotonically decreasing values of (i't), 
but the selection improves their average robustness significantly, above or up to the same level as the networks with 
(fc) = 2.5. The robustnesses of the initial state, 1^8(0)1 and the final state, I's*, (Figs, [^e) and (f)) follow the same 
trend. The robustnesses of the principal basin, v^, and state space, vg,, (Figs. [6l[g) and (h)) do not show as much 
improvement with selection. This implies that the effect of selection is more local (on the specific gene-expression 
trajectory) than global in terms of the change in dynamics. Our results above agree with the findings by Ciliberti 
et al. [l3|, who used very similar parameters to measure robustness to noise. 

We note that most of the robustness measures discussed above do not increase monotonically with (fc) for the viable 
networks (i.e., before selection). For example, (i?^) for viable networks has a maximum around (/c) ~ 5 (Fig. [6|a)). 
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FIG. 6: Mean mutational robustness {R^l), transient time (r), normalized size of the principal basin (B), robustnesses of the 
trajectory {vt), initial state (!/s(o)), final state (fs*), principal basin (^^b), and state space (fs) before (viable) and after evolution 
(evolved) for A^ = 10 and (k) = 1, 2.5, 5, 7.5 and 10 averaged over 10,000 networks. The error bars are equal to one standard 
deviation (not standard error). B is normalized by fl/2 — 2 ^^ . The lines connecting the symbols are guides to the eye. 

This probably indicates that these quantities depend on more than one parameter (Hke the principal basin size and 
the magnitude of the damage spreading), some increasing and some decreasing with increasing (fc). For the evolved 
networks, however, (i?^), (B), and (i^t) increase monotonically with (fc), at least for the range tested in this paper. 
This indicates that networks of higher connectivity are more evolvable when mutational robustness is considered. 
However, these networks are more chaotic on average compared to the ones of low connectivity, regardless of selection. 

D. Correlations between Change in Mutational Robustness and Other Parameters 

The recent theoretical and computational work on gene regulatory networks [10, UM indicates that there is a strong 
link between mutational robustness and robustness to noise. In order to see the relation between these two quantities, 
we calculated the correlation coefficients between the change in mutational robustness (Ai?^) and changes in different 
measures of robustness to noise, as well as transient time and the damage spreading as functions of linkage density. 

Changes in the transient time and mutational robustness are highly correlated in agreement with the earlier studies 
[43|. The principal basin size does not seem to have much effect on i?^, as the correlation between their changes is 
quite weak (Fig. JZUb)). Although a larger principal basin size means greater robustness to noise, this result suggests 
that it is not strongly selected for. Changes in robustnesses of the initial state, final state, principal basin, and state 
space are also correlated with AR^ (Figs. jTlJc), (d), (e), and (f)). The weak correlation between the changes in the 
damage spreading and mutational robustness (Fig. |7l[h)) imply that it is not the change in dynamical behavior that 
brings mutational robustness. 



DISCUSSION 



In this paper, we have analyzed changes in state-space properties of model genetic regulatory networks under 
selection for an optimal phenotype. Both numerical stability analysis and the state-space statistics indicate that the 
difference between the global dynamical properties of mutationally robust networks that have undergone selection and 
their random ancestors are quite small. Furthermore, the correlation between the changes in the damage spreading 
and the mutational robustness is weak. Therefore, changes in the global dynamical properties do not seem to be 
responsible for the increase in mutational robustness after selection. 

Dynamics of many random threshold networks, as well as Random Boolean Networks, depend largely on their 
connectivity distributions. Variants of the Random Threshold Networks used in this paper have been shown to have 



10 




FIG. 7: Correlations between change in state-space parameters and change in robustness, calculated using 10,000 individual 
evolved networks with A'' = 10 and (A;) = 1, 2.5, 5, 7.5, and 10. The P-value is much smaller than 0.01 for all correlation 
coefficients with an absolute value above 0.03. Ar, AS, Az^t, Az^s(o), Az/s*, Ai-b, Afs, and Adt+i denote changes in normalized 
size of the transient time, size of the principal basin, robustnesses of the trajectory, initial state, final state, principal basin, and 
state space, and damage spreading, respectively. The dynamical robustness of the trajectory includes the contribution from 
the final state. The lines connecting the symbols are guides to the eye. See text for details. 



a 2, depending on the model details [26 
,\EM just as RBNs SQIll. 



|34| - Essentially, RTNs have a chaotic phase for 



a chaotic phase above 
sufficiently large (fe) [2£ 

The connectivity of a network affects the evolvability of its mutational robustness, as well as its dynamical character. 
For viable (essentially random) networks, the mutational robustness is very similar for all connectivities [43|. For the 
evolved networks, however, it increases monotonically with increasing connectivity, creating drastic differences for 
large (fc). The dynamical robustness of the gene-expression trajectory, the initial state, and the final state follow 
similar trends. These results clearly indicate that the stability of the system as measured by the damage spreading 
does not capture its dynamical characteristics in this context. 

As pointed out in a recent paper by Ciliberti et al. [10||, selection decreases the transient time by pic king the 
"proper" interaction constants to construct a shorter (or direct) path from s(0) to s*. Both Ciliberti et al. [10| and 
Kaneko [2l| stated that mutationally robust networks are the ones that have found a path for the gene expression 
trajectory at a safe distance from the basin boundary, so that small perturbations cannot kick them into a different 
basin. Thus, even if the selection operates only on the stability of the stationary gene-expression pattern, robustness 
to mutations intrinsically requires stability of the gene-expression trajectory against small perturbations. There is also 
some experimental evidence supporting the association between genetic and non-genetic change [33,[43|- Although it 
might seem like robustness to noise evolves as a by-product of robustness to mutations, the converse case is also true 
21 1 . Thus, robustness to noise and robustness to mutations seem to evolve mutually when certain conditions are met 

m 



The analysis of the cell-cycle regulatory network of budding yeast provides further evidence for the chaotic behavior 
of gene networks |28l. l29l|. The simplified form of this network has 11 nodes (genes or proteins), one checkpoint (an 
external input, in this case, the cell size), and 34 links including self-degrading interactions. Using a dynamical model 
similar to the one used in this paper, Li et al. [23| showed that the stationary state of the network has a basin occupying 
86% of the state space. Lau et al. [23| studied an ensemble of networks that can perform the same function, i.e., 
the 12-step sequence of transitions in the expression trajectory, and found that these functional networks have larger 
basins for the stationary state (consequently, broadly distributed basin sizes), fewer attractors, longer transient times, 
and a larger damage-spreading rate compared to their randomized counterparts. They concluded that those dynamical 
features emerge due to the functional constraints on the network. Here, we showed that those features, which are 
signs of chaotic dynamics, can arise under the presence of structural perturbations (mutations) if the connectivity of 
the network is large enough, even when the constraints on the function are minimal, i.e., when the only selection is 
on the phenotype. Although the length of the gene-expression trajectory of the yeast cell-cycle network needs further 
explanation in terms of mutational robustness, it appears like chaotic dynamics may be a design principle underlying 
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seemingly "boring" and ordered behavior generally seen in models of gene regulatory networks, where a simple cascade 
of expression terminates at a stationary state [40| • 

The effect of network topology on evolvability of robustness is another aspect of the problem, which we do not 
discuss in this paper [19l . l20|, |33|. However, we would like to point out that recent studies indicate that certain 
topological features, such as connectivity, are not very crucial in determining the response of cellular networks to 
genetic or non- genetic change [la, [33|, and there may be other factors shaping their topological structure |6l.[l5l|. 

Our results also imply that the "life at the edge of chaos" hypothesis [231 , which suggests adaptability (evolvability) 
is maximized in the critical regime does not seem to be necessary, at least not to explain evolvability of robustness to 
mutations and noise. Indeed, recent studies concerning dynamics of genetic regulatory networks do not indicate any 
special feature brought by criticality [2| • The "edge of chaos" concept was primarily developed to describe the phase 
in which cellular automata can perform universal computation |27l. l32l|. and it may not be related to dynamics and 
evolution of biochemical pathways as was once thought. 

To summarize, our study indicates that conventional measures of stability may not be very informative about 
robustness to mutations or noise in gene regulatory networks when one considers the steady gene-expression pattern 
as the robust feature of the network. Also, the dynamics underlying the simple gene-expression trajectories can be 
very rich, reflecting a complex state-space structure. 

Supplementary material 

The online version of this article contains supplementary material: Linearly binned histograms of the data rep- 
resented in Figs. [4] and [H probability densities for the Hamming distance between the initial and final states, the 
magnitude of the damage spreading before and after evolution, and probability densities for the latter. 
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